
library(haven)

experiment <- read_dta("experiment_replication_data.dta")

names(experiment)

#voting

experiment$akp[experiment$vote == 1] <- 1
experiment$chp[experiment$vote == 2] <- 1

experiment$akp[is.na(experiment$akp)] <- 0
experiment$chp[is.na(experiment$chp)] <- 0


#generate experiment variables

experiment$ott <- NA
experiment$ott[experiment$treatment==1] <- 1
experiment$ott[is.na(experiment$ott)] <- 0

experiment$kem <- NA
experiment$kem[experiment$treatment==2] <- 1
experiment$kem[is.na(experiment$kem)] <- 0

experiment$con <- NA
experiment$con[experiment$treatment==0] <- 1
experiment$con[is.na(experiment$con)] <- 0


experiment$ts[experiment$treatment == 0] = "Control"
experiment$ts[experiment$treatment == 1] = "Ottoman"
experiment$ts[experiment$treatment == 2] = "Kemalist"

experiment$ts = factor(experiment$ts, levels = c("Control", "Ottoman", "Kemalist"))

library(ggplot2)
library(ggrepel)
library(Rmisc) 
library(psych)
library(gridExtra)
library(stargazer)
library(jtools)
library(ggstance)
library(huxtable)

#sample characteristics (Table 3)

names(experiment)

describe(experiment[,c(13, 12, 14, 15, 16)])
describe(subset(experiment[,c(13, 12, 14, 15, 16)], experiment$ott == 1))
describe(subset(experiment[,c(13, 12, 14, 15, 16)], experiment$kem == 1))
describe(subset(experiment[,c(13, 12, 14, 15, 16)], experiment$con == 1))

#multinomial logistic regression

library(nnet)

test <- multinom(treatment ~ age + female + education + akp + chp, data = experiment)

z <- summary(test)$coefficients/summary(test)$standard.errors
z

p <- (1 - pnorm(abs(z), 0, 1)) * 2
p

library(lmtest)

lrtest(test)

summary(test)

#experimental analysis

#scale building

experiment$popr <- experiment$pop2  + experiment$pop3 + experiment$pop4  +  experiment$pop5 + experiment$pop8

experiment$popr <- (experiment$popr+10)*5

summary(experiment$popr)

names(experiment)

##Table 4

cor(experiment[,c(5, 6, 7, 8, 11)], use = "pairwise.complete.obs")

describe(experiment[,c(5, 6, 7, 8, 11, 21)])



#Mean plot (Figure 1)


populism <- summarySE(experiment, measurevar="popr", groupvars="ts", na.rm = T)

ggplot(populism, aes(x=as.factor(ts), y=popr)) +
  geom_errorbar(width=.1, aes(ymin=popr-ci, ymax=popr+ci), size=0.8) +
  geom_point(shape=21, size=3, fill="white") +
  scale_y_continuous(breaks = seq(60,80,by=5), limits = c(60,80)) +
  #geom_hline(yintercept = 70.69)+
  labs(x="Treatment groups") +
  labs(y="Populism") + theme_bw() +
  theme(axis.text.x = element_text(face = "bold")) +
  theme(axis.title.y = element_text(face ="bold")) +
  theme(axis.text.y = element_text(face = "bold")) +
  theme(axis.title.x = element_text(face ="bold")) +
  theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor = element_blank()) +
  theme(axis.text.x=element_text(colour="black")) +
  theme(axis.text.y=element_text(colour="black"))
#+ theme(text = element_text(size = 15))

library(ggtern)

ggsave("populism.tiff", dpi = 1200, width = 625, height = 365, units = "mm", scale = 0.3)

#ggsave("populism1.tiff", dpi = 1200, width = 625, height = 365, units = "mm", scale = 0.3)

#anova test

popanova <- aov(popr ~ ts, data = experiment)

summary(popanova)

TukeyHSD(popanova)

#regression

model1 <- lm(popr ~ ott + kem, data=experiment)

summary(model1)


model2 = lm(popr ~ ott*akp + kem*akp, data = experiment)

summary(model2)

model3 = lm(popr ~ ott*chp + kem*chp, data = experiment)

summary(model3)


export_summs(model1, model2, model3, scale = F, ci_level = 0.90, to.file = "docx", file.name = "populism.docx")

#Average treatment effects (Figure 2)


confint(model1)

coef = c(3.664, 2.963)

cilow = c(0.261, -0.341)

cihigh = c(7.066, 6.266)

variable = c("Ottoman", "Kemalist")

label = c("3.664*", "2.963+")

data = data.frame(variable, coef, cilow, cihigh, label)

ggplot(data, aes(x=as.factor(variable), y=coef)) +
  geom_errorbar(width=0, aes(ymin=cilow, ymax=cihigh), size=0.8) +
  geom_point(shape=21, size=3, fill="white") +
  geom_hline(yintercept = 0)+
  labs(x="") +
  labs(y="Coefficients") + 
  geom_text(aes(label = label), hjust = -0.15, nudge_x = 0.1, size = 3.5) +
  theme_bw() +
  theme(axis.text.x = element_text(face = "bold")) +
  theme(axis.title.y = element_text(face ="bold")) +
  theme(axis.text.y = element_text(face = "bold")) +
  theme(axis.title.x = element_text(face ="bold")) +
  theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor = element_blank()) +
  theme(axis.text.x=element_text(colour="black")) +
  theme(axis.text.y=element_text(colour="black")) +
  theme(text = element_text(size = 12)) +
  coord_flip() +
  scale_y_continuous(breaks = seq(-2,8,by=2), limits = c(-2,8))

ggsave("regression.tiff", dpi = 1200, width = 625, height = 365, units = "mm", scale = 0.3)
